Grand-canonical variational approach for the t — J model 
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Gutzwiller-projected BCS wave function or the resonating-valence-bond (RVB) state in the 2D 
extended t — J model is investigated by using the variational Monte Carlo technique. We show that 
the results of ground-state energy and excitation spectra calculated in the grand-canonical scheme 
allowing particle number to fluctuate are essentially the same as previous results obtained by flxing 
the number of particle in the canonical scheme if the grand thermodynamic potential is used for 
minimization. To account for the effect of Gutzwiller projection, a fugacity factor proposed by 
Laughlin and Anderson few years ago has to be inserted into the coherence factor of the BCS state. 
Chemical potential, particle number fluctuation, and phase fluctuation of the RVB state, difficult 
or even impossible to be calculated in the canonical ensemble, have been directly measured in the 
grand-canonical picture. We find that except for La — 214 materials, the doping dependence of 
chemical potential is consistent with experimental findings on several cuprates. Similar to what has 
been reported by scanning tunneling spectroscopy experiments, the tunneling asymmetry becomes 
much stronger as doping decreases. We found a very large enhancement of phase fluctuation in the 
underdoped regime. 
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PACS numbers: 71.27+a,71.10- 



I. INTRODUCTION 



Since the pioneering work done by Andersoni, the su- 
perconducting (SC) state of high- Tc cuprates has been 
successfully described by the so-called d-wave resonating- 
valence-bond (d-RVB) wave function with a fixed number 
of particles, 



K-RVB)=PNpG\^d-BCs) 



(1) 



where the Gutzwiller projection operator Pq restricts 
Hilbert space without doubly occupancy at each site. 
Pn^ is the projection operator onto the subspace with 
Nf. electrons. According to Eq.(IT]), the SC wave function 
with variable particle numbers seems to be naturally ob- 
tained by taking away Pm^ ■ However, recently Anderson 
has emphasized that a fugacity factor should be inserted 
in front of the coherence factor Uk in Eq.(IT]) besides re- 
moving Pn^^'^- This is the same idea as what Laughlin 
proposed for the Gossamer superconductivity-.. By using 
Gutzwiller approximation (GA) Edegger et al^ have also 
discussed the necessity of introducing the fugacity factor 
in the grand canonical wave function. But there is no ex- 
act numerical result to veriiy it. Whether GA provides 
an accurate fugacity factor is still an open question. 

Although Eq.dl]) has been used to explain several im- 
portant features of high- Tc cuprates like the ground- 
state phase diagram^"^9., interplay between antiferromag- 
netism and superconductivitjJ^^— , existence of stripe 
statesi^rJS, and anomalous spectral weights of low-lying 
excited states^"—, etc., only few numerical works on 
the grand-canonical c?-RVB state have been reported so 
far. As far as we know, Yokoyama and Shiba is the 
first to have carried out the variational Monte-Carlo 
(VMC) calculations with non-conserving particle num- 
bers in strongly correlated Hubbard model almost two 



decades ago2^. By using particle-hole transformation, 
the calculation can be efficiently performed. However, 
they did not consider the fugacity factor in the calcula- 
tion. They just briefly mentioned how to introduce an 
additional variational parameter a in front of the coher- 
ence factor Wfc to control the average particle number. 
We also notice a will become very large near half-filling 
and that causes serious numerical difficulties. Thus, it is 
important to re-examine this approach with the fugacity 
factor added in front of Uk instead. 

Using this grand-canonical wave function we can ex- 
amine several important physical quantities that were not 
able to obtain by fixing number of particles. For instance, 
Anderson and Ong'^ showed the importance of the fugac- 
ity factor by using the GA to calculate the famous asym- 
metric tunneling conductance observed by scanning tun- 
neling spectroscopy (STS)^^"^^. Although in our earlier 
studies^ using Eq.(IT]), we have shown numerically the 
asymmetry arise from strong correlation effect, this con- 
clusion needs to be verified in the grand-canonical calcu- 
lation. Since the conductance involves the ground-state 
energy, excitation spectra, and spectral weights, now the 
excitation involves Bogoliubov quasi-particles unlike the 
excitation of Eq. ([Ij which has quasi-particles with a def- 
inite charge. In addition, the SC order can now be cal- 
culated explicitly instead of using the long-range pair- 
pair correlation function to get an estimate. The par- 
ticle number fiuctuation or the phase fluctuation in the 
strongly correlated SC state can be calculated directly 
also. Furthermore, the doping dependence of chemical 
potential which was reported by experiments^^, can be 
determined directly in this grand-canonical scheme. 

Below we shall flrst focus on computing several phys- 
ical quantities using the d-RVB state with fluctuating 
particle numbers, such as the nearest-neighbor SC order 
parameter Ad(= X]i(^!t^l+S4.))' particle number flue- 
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tuation AA''e(= a/ {N"^) — (Ne)^), chemical potential ^g, 
and spectral weight Z^^, etc. First of all, the validity 
of d-RVB wave function including a fugacity factor in 
the grand-canonical picture is numerically confirmed, as 
we find the equivalence for the ground- and excited-state 
features between the canonical and the grand-canonical 
ensemble by optimizing the grand thermodynamic po- 
tential F instead of internal energy E. Next, the doping 
trends of chemical potential and tunneling conductance 
are evaluated, which is in agreement with the experimen- 
tal observations. Both SC order parameter and particle 
number fluctuation show the similar doping dependence 
with a dome-like shape observed in the SC phase of high- 
Tc phase diagram. Finally, we shall examine the effect of 
strong correlation on the phase field of the SC order pa- 
rameter. We find that in the underdoped region, the 
strong correlation effect greatly enhances phase fluctua- 
tion in contrast with the prediction of BCS theory. 



II. VARIATIONAL MONTE CARLO METHOD 

The Hamiltonian for the extended t — J model on a 
two-dimensional square lattice is given by 



H 



h.c. 



{i,3) 
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(2) 



where the hopping terms — t, t', and t" for sites i 
and j being the nearest-, the second-nearest, and the 
third-nearest-neighbors, respectively. Other notations 
are standard. The SC wave function is of the form. 



kf^'-kj. 



|0), (3) 



where the coefficients Wk and Uk are not necessarily the 
coherence factors in the BCS wave function. 

To analyze further we turn to introduce the framework 
of VMC approach. The expectation value of the Hamil- 
tonian H using the d-RVB state is evaluated as^ 



E EE 



d-RVB 
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d-RVB 



(4) 



Here p/j represents the Monte-Carlo sampling weight for 
the configuration j3 defined as 



I* 



d-RVB l\ — 



det 



(5) 



In the canonical ensemble, the matrix elements 
is given by 



(6) 



where Ri-|- and Hji is the position vector of the i-th up- 
spin electron and the j-th down-spin electron in the con- 
figuration /3, respectively. In Eq.([3]), the particle number 
fluctuation results in the unknown matrix size bringing 
technical difficulty to VMC calculations. To resolve this 
problem, a spin-dependent particle-hole transformation 
has to be introduced. 

Following Yokoyama and Shiba^^ we introduce the par- 
tial particle-hole transformation to change the original 
representation (c) to a new representation (df ) expressed 
as 



Cif fi, 

Cii -> d], 



(7) 



where only the down-spin electrons are transformed. 
Here we introduce two different particles, d and /, instead 
of down- and up-spin electrons. Both operators fi and di 
also satisfy the anti-communication relation. Thus, three 
possible Fock states at each site can be transformed in 
the following way. 



|0)(c) 
I i)(c) 
lt)(c) 



midf) 

|0)(./) 
\df)(df)- 



(8) 



The subscripts indicate different representations. Now 
Eq.dSl) can be transformed into the representation (df). 



PGn("k + «k4/_ki) I0)(c) 
k 

k 



(9) 



The Gutzwiller projection operator Pq in the repre- 
sentation (df) restricts the Hilbert space to the three 
possible states shown in Eq.®. Eq.® displays a quan- 
tum state that the total number of d and / particles is 
fixed to the lattice size N. Thus there is no total parti- 
cle number fluctuation in the representation (df). This 
conservation can be understood from Eq.® as well. It 
suggests that for total Sz = the number of empty and 
doubly-occupied sites are always equal in the represen- 
tation (df), implying the total number of d and / par- 
ticles is exactly equal to N. However, the number of d 
and / particles can vary even though the sum is fixed. 
This fluctuation replaces the particle number fluctuation 
in the original representation (c). The coherence factors 
in Eq.(|9l) determine the number of d or /. Notice that 
in the representation (df) the average number difference 
between the particles d and / is equal to doping density 
S. 

In the canonical ensemble we usually have two Monte 
Carlo processes: hopping and exchanging. In other 
words, an up or down spin can hop to a hole site and two 
anti-parallel spins can exchange with each other. We can 
generate all states by applying one or both of the pro- 
cesses sequentially. To connect the Hilbert spaces with 
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different particle numbers in the grand-canonical scheme, 
however, we have to consider a new process to create or 
annihilate pairs. To change from a configuration with 
Ne particles to another with ±2 particles, we ran- 
domly choose two different sites i and j. If the sites i 
and j have opposite spins, we destroy the pair at these 
two sites. Conversely, if both sites i and j are empty, 
we create a pair. In the representation (df) these pro- 
cesses also can be easily implemented by changing two 
sites both with a single d particle to an empty site and a 
site doubly occupied with d and / particles or vice versa. 



III. THE d-RVB WAVE FUNCTION IN THE 
GRAND-CANONICAL ENSEMBLE 

According to the earlier results from GA^, they showed 
that Gutzwiller projection not only imposes a local con- 
straint on the electron number at each site, but glob- 
ally influences the total number of electrons. In the 
following, we will provide another argument to obtain 
the d-RVB wave function with non-conserving particle 
numbers. Since there exists particle number fluctuation 
A7Ve(c>c Vn) stemming from the SC order, the electron 
number in the SC ground state will have a distribution, 
p{Ne). The distribution function for the wave function 
without Gutzwiller projection, po{Ne), has an approx- 



imate Gaussian form e 



')V(AAr(' 



centered at 



the most probable number of electrons ivi°'' with a width 
ATVi"'. In the thermodynamic limit, Ne*^^ is the aver- 
age number of electrons. We define = N {1 — S) and 
ANe^^ = aoVW so that the variable Nf, in po{Ne) can be to 
changed into the doping S, 



(10) 



The distribution function of the Gutzwiller-projected 
wave function p{S), reduced by a weighting factor 
Ws, deviates from po{S) due to the constraint of no- 
doubly occupancy^. This weighting factor can be es- 
timated roughly by counting the number of configura- 
tions in the phase space. Before Gutzwiller projection, 
there are Cf^(^i_s)/2 choices for the up-spin or down- 
spin electrons to give us a total configuration number 
Nb (— {C^,-i_;i)/2)'^^- After the projection, the total 

-,N(l+6)/2 



number is changed into Na [= C^(i_5)/2 ' '^n{i-S)/2 
Thus, the weighting factor can be written as 



[iN{l + S)/2)\Y 

m (Nsy. 

[{N{l + 6)/2)lY 
[{N/2)lf {NSy. 



Wn 



(11) 



Now we assume N is very large. By using Stirling's 
approximation (TV! w N^/e^), Eq. pi|) can be reduced 
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FIG. 1: (Color online) (a) The average doping density and 
(b) its fluctuation as a function of the variational parameter 
A„ for the wave function of Eq.® using the BCS coherence 
factors. Results for a cluster of 10 x 10 and 18 x 18 are shown 
in black and red, respectively. Other parameters are set to 
be pv = t'^ = t'v = 0. The blue line in (b) is the d-wave BCS 
result. 
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Eas. (fTO| and (fT2|) lead to the distribution function p{5) 
p{d) = poiS) ■ Ws 
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This function is very different from po{S). Not only the 
maximum of the distribution function p{S) is shifted from 
So to the higher doping density 6, the width also becomes 
considerably narrower. In the thermodynamic limit, S in 
the distribution function is just the average doping den- 
sity {6). Thus, the doping density in the grand-canonical 
ensemble is determined not only by variational parame- 
ters but also the Gutzwiller projection operator. Similar 
conclusions have also been discussed previously^^. 

We can further study Ea. p3)) by carrying out the nu- 
merical calculation using the VMC method. We consider 
the simple BCS wave function by replacing {tk and Vk in 
Eq.® with Uk and Uk, respectively. Wk and Vk are the 
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BCS coherence factors, defined by 



\ 



Vk = Sgn{Ak) 



\ 



(14) 



where 



Efc = ^2{coskx + cosky) ^ At'^coskxCosky 



— (cos2kx + cos2ky) 
2Ay {coskx — cosky) . 



(15) 



Here At,, i^, and t" are variational parameters. For 
illustration we will consider the half-filling wave function 
with the variational parameters: = t'^ = t'l — Q. The 
only parameter left is Au which varies from to 1. In 
FigUJa) and (b), we show the average doping density 
((5)(= 1 — Ne/N) and the particle number fiuctuation 
ANe/VN as a function of A„, respectively. The result of 
the d-wave BCS state without projection is also shown in 
FiglUb). Clearly as discussed above the introduction of 
projection has greatly changed the distribution function 
of particle number with a larger average doping density 
and smaller fiuctuation. 

Now we can construct the d-RVB wave function in 
the grand-canonical ensemble by taking into account this 
change of distribution function. Eq. ([T2|) suggests that the 
important doping dependence of the distribution function 
in the presence of the Gutzwiller projection operator Pq 
is the factor g~^'^, where g = y^-, = N — N^, and 
Ne — J2i a ^icr^i<^- The Operator g^^'^ mimicking the 
effect of Pg tends to increase the doping density away 
from half-filling. To balance this effect, we have to place 

the operator y^^'' in front of the d-RVB wave function. 
Thus, the d-RVB state in the grand-canonical scheme can 
be written as 



I*: 



-RVBI 



+ t^k4tclki) |0)(c) 



Pg\^1-bcs) 



(16) 



Now g can be seen as a variational parameter like /j,^ , , 
t", and At,. Eq. (fT6l) with the fugacity factor is exactly 
the wave function proposed by Anderson^ and Laughliiii. 
Previously the fugacity factor g was estimated to be 
by using GA. 

To obtain the ground state in the grand-canonical case, 
we have to optimize the grand thermodynamic potential 
F — E ~ figNe instead of the internal energy E with 
a fixed chemical potential /ig. Here is the average 
number of electrons for each ^g. Notice that fig is differ- 
ent from the variational parameter /z„ in the c?-RVB wave 



■ • — g (with Jastrow) 




FIG. 2: (Color online) The doping dependence of (a) the vari- 
ational parameter g and (b) the optimized energy per site for 
a 12 X 12 lattice. In (a). Solid (Empty) circles indicate the fu- 
gacity factors g with (without) Jastrow factors. The blue line 
shows the renormalized Gutzwiller's factor ^f^- The green 
lines are the guide to the eyes. In (b), black squares (Red cir- 
cles) indicate the d-RVB wave function with Jastrow factors 
in the (grand-) canonical ensemble. Error bars represent the 
average density. Inset: Solid (Empty) squares represent the 
d-RVB wave function with (without) Jastrow factors in the 
canonical ensemble. The bare parameters in the Hamiltonian 
are set to be {t',t",J)/t = (-0.3,0.15,0.3). 



function. To have a lower ground-state energy, for all the 
numerical results discussed below we consider a modified 
d-RVB wave function Pjl'^d.-Rvs) 'where we introduce a 
hole-hole repulsive Jastrow factor P j^^i'^^" — : 



P.I 



i<j 



with 



(17) 



(18) 



Here n f = 1 — c|g.Cio-- The three parameters vp with 
/? to be the nearest, second nearest, and third nearest 
neighbors are for short-range hole-hole repulsion if these 
values are less than 1. The factor rg is for long-range 
correlations and it is repulsive if a is positive. L is 
the linear scale of the lattice. In Figl^Ja), we obtain 
the doping dependence of the fugacity factor g from the 
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grand-thermodynamic-potential optimization. Interest- 
ingly, since Jastrow factors will reduce the probability 
that holes come closer, the fugacity factor of the modi- 
fied d-RVB state is greatly enhanced as increasing dop- 
ing. If we do not consider Jastrow factors, the doping de- 
pendence of g is similar to the renormalized Gutzwiller's 
factor 2(5/(1 -f 5). Therefore, although the GA result is 
not exact, it is a reasonable approximation. 



IV. QUASI-PARTICLE ENERGY DISPERSION 
AND SPECTRAL WEIGHT 

In Figl^Ib), the optimized internal energy per site is 
plotted as a function of the doping density for both 
canonical (black squares) and grand canonical (red cir- 
cles) cases. The results are essentially indistinguishable 
for a 12 X 12 lattice. Thus the ground-state phase diagram 
obtained in the canonical ensembl o"'^^'^^ will be essentially 
the same as the grand-canonical scheme. This result is 
expected as the calculation of internal energy only in- 
volves states with the same number of particles^. Ad- 
ditionally, the energies obtained with the Jastrow factor 
are quite lower than without the factors for the canonical 
case, as shown in the inset of Figl^Kb). The fugacity fac- 
tor g of the modified d-RVB state is very different from 
the case without Jastrow factors shown in Figl^^a). 

It should be noted that we can also use the results 
in Figl2Kb) obtained by Eq.([T]) to calculate the relation 
between chemical potential and average number of parti- 
cles as /ig = dE/dNe- Completely consistent results are 
obtained. The excellent agreement between the energies 
calculated by the grand-canonical and canonical schemes 
is the most important numerical result of this paper to 
firmly establish the validity of inclusion of the fugacity 
factor g in Eq. (flB]) and our Monte Carlo algorithm. Now 
we can start to calculate excitation spectra and other 
physical quantities in the grand-canonical ensemble. 

Not only we shall consider the energy dispersion 
of the excitations but also the spectral weight be- 
low. These results could be compared with the mea- 
sured STS. The simplest way to construct a single- 
particle excitation from the d-RVB state is to bring the 



quasi-particle creation operator 7j,g.( 
■■c-k-c) into play. 



(19) 



Excitation energies cannot be calculated from the in- 
ternal energy difference — Eq as in the canonical 
ensemble^^^, where Eq is the ground-state energy. Here 
the chemical potential is fixed, hence we must calculate 
the grand thermodynamic potential for the ground state 
and excited states, _Fo and Fk, respectively. 



0.4-1 
0.3 - 

O0-2- 
0.1 - 
0.0 



(a) 



A 



Canonical 
► — Grand-Canonical 



0.3 -I 



> 



.0.2 - 



0.1 - 



0.0 



(b) 




<5>=0.I10 
<5>=0.077 
<5>=0.039 



-1.0 



-0.5 



0.0 

v/t 



0.5 



1.0 



FIG. 3: (Color online) G{V) for the d-RVB state versus the 
bias V. (a) Canonical (Grand-canonical) results for hole 
doping 5 (chemical potential pg) fixed to 0.125 (1.68f) in 
12 X 12 lattice, (b) Grand canonical results for three chem- 
ical potentials ^g, 1.68t (black-squares), 1.9t (red-circles), 
and 2. It (green-triangles). V is negative (positive) for re- 
moving (adding) one electron in 20 x 20 lattice. The bare 
parameters in the Hamiltonian are set to be {t' ,t" , J) /t = 
(-0.3,0.15,0.3). 



It is noticed that due to the BCS coherence factors Uk 
and Vk, the excited state l^'f^,.) will have fewer (more) 
average particle numbers below (above) the Fermi level 
than the ground state l^'^-itys) i^^^ shown). Here we 
have made an important assumption that these partic- 
ular states Eq. p^ have little overlap with other states 
with same momentum and spin. This is probably valid 
for low-energy excitations less than the gap energy. 

The above quasi-particle states are then used to calcu- 
late the spectral weight for adding (removing) one elec- 
tron with momentum k {—k) and spin a (— cr) to the 
ground state as defined 
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Similar to what we have done for the canonical case^^, 
we could also calculate the tunneling conductance for 
the grand-canonical wave function with the bias given by 
V = Fk — Fq. For the numerical calculations, we define 
the tunneling conductance as 
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(20) 
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FIG. 4: (Color online) (a) The doping dependence of chem- 
ical potential fig in a 12 x 12 lattice. Red circles represent 
the d-RVB state with Jastrow factors in the grand-canonical 
ensemble. Error bars indicate ANe/N at each average dop- 
ing, (b) The doping dependence of chemical-potential shift 
from experiments in several cuprate materialsSS, compared 
with the region shown by the green-dashed frame in (a). Red 
empty circles in (b), shifted along the vertical direction for 
comparison, are the data from the green-dashed frame in 
(a). The bare parameters in the Hamiltonian are set to be 
[t',t",J)/t = (-0.3,0.15,0.3). Here t is set to be 0.461^. 



0.10 n 0.3 n 




I ^ — . — 1 — . — 1 — . — 1 — . — 1 — . u.u ^ — i — I — i — I — i — I — i — I— 
0.0 0.1 0.2 0.3 0.4 0.0 0.1 0.2 0.3 0.4 



S 8 

FIG. 5: The doping dependence of (a) the SC order parameter 
Ad and (b) the scaled particle number fluctuation ANe/V^ 
using the d-RVB state in the 12 x 12 lattice in the grand- 
canonical ensemble. Error bars indicate ANe/N at each av- 
erage doping. The bare parameters in the Hamiltonian are 
set to be {t',t",J)/t = (-0.2,0.1,0.3). 



where AF = 0.28t (0.2t) for N = 144 (400) chosen is 
to reduce the effect due to the finite lattice size and a 
either up or down spin. Similar to Figl2jb), FiglSja) 
also shows except for high voltage the tunneling con- 
ductances are almost identical in the canonical and the 
grand-canonical ensemble. Once again, the equivalence 
between the canonical and the grand-canonical ensemble 
for the excitation spectra within the gap convinces us of 
the conclusion that the modified d-RVB wave function 
with the fugacity factor is the precise representation for 
the RVB-type state in the grand-canonical case. 

In addition, to numerically examine the particle-hole 
asymmetry of tunneling spectra in the grand-canonical 
ensemble, we present the doping dependence of G{V) in 
FigEJb). The asymmetry is clearly observed for all three 
doping densities. At <5 = 0.039, the total spectral weight 
for removing one electron, defined as the sum of Zf^^ 
over the Brillouin zone, is found to be about three times 
as large as the one for adding one electron. The gap 
value deduced from the width between peak positions 
decreases with doping. Note that the gap size is usually 
overestimated in the extended t — J models. However, 
the weight value obtained from the peak height enhances 
as increasing doping, apparently anticorrelating with the 
gap size. All of these features have been shown in the 
canonical ensemble^ and qualitatively consistent with 
the STS measurements^^. Although the details of the 
tunneling spectra in Ref. is not completely identical to 
that shown in Fig[3l^b), their important characteristics 
are similar. 



V. CHEMICAL POTENTIAL 



The doping dependence of fig calculated in the grand- 
canonical ensemble is plotted in FigUJ^a). It decreases 
monotonically with doping as expected from Fig|4ljb). 
Near half-filling, fig seems to increase greatly. The com- 
pressibility becomes extremely small as electron density 
approaches half filling, this is a consequence of the open- 
ing of a charge gap of a Mott insulator. This behavior 
was also observed in the calculation of Yokoyama and 
Shiba^^. But the quantitative behavior is different as the 
fugacity factor was not included in their grand-canonical 
wave function. We have also examined the doping de- 
pendence of fig for several bare parameters t' /t. The 
slopes for all t' /t are similar except for very large doping 
(not shown). Here, only the relative value of fig is mean- 
ingful. Thus we can shift fig to compare with experi- 
ments. FigureSl^b) shows chemical-potential shift for sev- 
eral cuprates observed by photoemission experiments—. 
Except for La — 214 samples the slope of chemical poten- 
tial as a function of S obtained by our calculation (see the 
green-dashed frame in FigHl^a)) is almost identical to the 
experimental data. This result is more quantitatively re- 
liable in comparison with the experiments than our previ- 
ous studies in the canonical ensemble^ in the low doping 
region. As for the inconsistency with La — 214 cuprates, 
a possible reason could be the existence of stripe order 
in those materials^^^. It will be interesting to inves- 
tigate chemical-potential shift by using the Gutzwiller- 
projected stripe wave functioi*i^ in the grand-canonical 
ensemble in the future. 



7 



VI. SC ORDER PARAMETER AND PARTICLE 
NUMBER FLUCTUATION 

In order to investigate the SC characteristics of the d- 
RVB wave function, we calculate the SC order parameter 
Ad and the particle number fluctuation AiVg impossi- 
ble to be derived in the canonical ensemble, as shown in 
FiglSJ^a) and (b), respectively. The SC order parameter 
Ad which is presumably proportional to the SC critical 
temperature Tc is plotted as a function of doping density 
in FiglSJa). The dome-like shape comparable to the SC 
phase in the cuprate phase diagrams is consistent with 
earlier VMC results in the canonical ensembleS^,. The 
density with the maximum value of order parameter is 
again much larger than the optimal doping density in 
cuprates. In the BCS theory the fluctuation of the par- 
ticle number, AN^, is proportional to the SC order pa- 
rameter. In FiglSjb), AN^/Vn is plotted as a function 
of doping density. Comparison between FiglSJa) and (b) 
shows that both AN^ and calculated by using the d- 
RVB state have the similar dome-like shape, but there is 
a significant difference at underdoped regime. This differ- 
ence is clearly due to the projection operator or the Mott 
physics as the particle number fluctuation is suppressed 
for low doping. 



VII. PHASE FLUCTUATION 

Before looking into the phase fluctuation A6 of the d- 
RVB wave function, we shall start with the c?-wave BCS 
state at flrst. To obtain useful information about the 
phase of wave functions, a "phase" operator is deflned 
by 



e 



A- At 
2iAo ' 



(23) 



where A^ = </5fecJ,^cLj,^ and the normalization Ag — 
|(At)|. Here we choose (pk — (Cfe-|-cLfe^)- Hence, any real 
wave function will lead to (0) = 0. With a little algebra, 
it is straight forward to write down the particle number 
fluctuation AN^ and the phase fluctuation AQ of the 
BCS state as 



AN, = 2 J2\uk\^\vk\\ 



AG = 



2l]fc 'Pk\uk\\vk\ 



(24) 



Then the uncertainty principle for particle number and 
phase can be easily derived by Cauchy-Schwarz inequal- 
ity: 



AA^eAe 



(J2kfk\uk\\vk\f 



> 1. 



(25) 




1/(NA0 ) (%)^o' 
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FIG. 6: The doping dependence of the particle number fluc- 
tuation AN^, the phase fluctuation AQ, and their product 
using (a) the d-wave BCS state and (b) the d-RVB state with 
Jastrow factors in the grand-canonical ensemble. The d-wave 
BCS state is obtained by solving the BCS Hamiltonian with 
bare parameters {t' ,t" ,V) /t = (—0.2,0.1,6) self-consistently. 
V is the pairing interaction strength. The blue dashed line 
denotes 1. (c) The doping dependence of AA^e (squares), AO 
(circles), and AA^^AO (triangles) using the d-RVB state with 
(filled) and without (empty) Jastrow factors, (d) The dop- 
ing dependence of ^^^^ with (filled diamonds) and without 
(empty diamonds) Jastrow factors. Black (Red) color indi- 
cates A*' = 144 (256). All the results are obtained with pa- 
rameters {t',t",J)/t = (-0.2,0.1,0.3) for a 12 X 12 lattice 
except specially mentioned results in (d). 



It can be proved that AA^g A8 in terms of the definition of 
Eq. ([23|) is exactly equal to 1 in BCS theory. The doping 
dependence of the particle number fluctuation and the 
phase fluctuation in the c?-wave BCS case are shown in 
FigUa). 

Although it is impossible to evaluate the phase fluctua- 
tion for the wave function of Eq. ([T]) with a fixed number 
of particles, it is straight forward for the c?-RVB wave 
function by calculating the following quantity 



Ae = 



e- e 



AtAtl^s 



-RVB/ 



2AI 



(26) 



In addition to Aq, there are two quantities to be calcu- 
lated by means of the VMC approach: one is the pairing 
correlation operator A^A and the other A^At. Using 
the representation (df), we can calculate both quantities 
directly. 

In Figinjb), we show the doping dependence of the 
particle number fiuctuation and the phase fiuctuation us- 
ing the d-RVB wave function. As mentioned before, the 
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Gutzwiller projection operator Pq suppresses the par- 
ticle nuniber fluctuation as shown in FiglSJ^b), especially 
for the underdoped region. On the other hand, the phase 
fluctuation is greatly enhanced for < (5 < 0.15, and a 
much weaker dependence for doping greater than 0.15. 
Although the fluctuation behavior is approaching BCS 
results in the overdoped region, it is still much stronger. 
This huge enhancement of A0 is clearly due to the strong 
correlation effect of the Gutzwiller projection operator. 
This result suggests that there may exist a strong phase- 
fluctuating state which is again consistent with experi- 
ments in the underdoped cuprate compounds^i^ and a 
theoretical analysis^. Another interesting result is the 
large enhancement of AiVeAO. Instead of having the 
value of 1 as the BCS state, it seems to approach infin- 
ity at very low doping. This is mainly due to the strong 
increase of phase fluctuation as doping decreases. 

FiglSlJc) shows the doping dependence of ANg 
(squares), A6 (circles), and AA^eA6 (triangles) using 
the d-RVB state with (filled) and without (empty) Jas- 
trow factors. As discussed before, the Jastrow factor 
suppresses the number fluctuation, hence it increases the 
phase fluctuation. But the substantial increase of the 
phase fluctuation at low doping region is quite surpris- 
ing since the energy difference between these two states 
is quite small as shown in the inset of Figl^Kb). Their 
product ANgAQ for the case with Jastrow factors also 
deviates farther from 1 near the underdoped regime. It 
indicates that the system in the underdoped region may 
have many low energy states with similar energy but dif- 
ferent properties. 

In FiglH^d), we find that the doping dependence of 
^Q^jv exhibits approximately a linear relation, shown by 
the empty diamonds, within 0.03 < S < 0.25 for the d- 
RVB state without including the Jastrow factor. Due to 
the finite size, it is difficult to get reliable results at ex- 
tremely low density but the results of two different clus- 
ter sizes are consistent. This result is in sharp contrast 
with BCS theory which has ^Ja^v proportional to the 
pairing gap instead of the doping density. After the Jas- 
trow factor is included, as shown by the filled diamonds 
in FiglHKd), phase fluctuation is enhanced. However, the 
linear dependence of doping density at low density still 
remains but with a smaller slope. 

Empirically the superfluid density of the hole-doped 
cuprates is small and proportional to the doping 
density^. For a SC state the phase stiffness is propor- 
tional the superfluid density. Although we do not have 
a direct proof of the relation between phase stiffness and 
^J^jv ; it is quite interesting that they both are propor- 
tional to the doping density. In addition, the slope for the 
wave function with the Jastrow factor is smaller than the 
one without the Jastrow factor, which may indicate that 
the SC ground state of the t — J model will have a very 
small superfluid density just as the hole-doped cuprates. 



VIII. CONCLUSIONS 



To summarize, we have studied the properties of the 
ground state and Bogoliubov quasi-particle states in the 
extended t— J model based on Gutzwiller-projected BCS 
wave function or d-RVB wave function in the grand- 
canonical ensemble. First of all, by using the phase space 
argument used in GA^, we have numerically demon- 
strated how to construct the correct d-RVB wave function 
with non-conserving particle numbers. A fugacity factor 
g in front of Uk in the d-RVB state is able to efficiently 
govern the distribution of empty sites. Our numerical 
calculations have shown the excellent agreement obtained 
for both the optimized grand thermodynamic potential 
and tunneling spectra G{V) between the grand-canonical 
and the canonical ensemble. It confirms the necessity 
and importance of the fugacity factor g in the d-RVB 
state in the grand-canonical ensemble as emphasized by 
Anderson^.. The enhanced tunneling asymmetry at low 
voltage in the underdoped region is again numerically 
demonstrated. 

In addition, as increasing doping, chemical potential 
fig calculated in the grand-canonical scheme monotoni- 
cally declines with the slope which is in good agreement 
with the experimental results"^". Almost zero charge sus- 
ceptibility near half-filling indicates the incompressible 
feature due to the Mott gap. Both the SC order param- 
eter and the particle number fluctuation have the dome- 
like doping dependence similar to the SC dome in high- Tc 
phase diagrams. The doping dependence of the SC order 
parameter is similar to that of the particle number fluctu- 
ation for large doping density as the BCS theory. But for 
low doping density, there is a significant difference due 
to the Gutzwiller projection operator. Furthermore, now 
we are able to directly calculate the phase fiuctuation in 
the grand-canonical ensemble. The Gutzwiller-projected 
wave function shows not only the smaller particle num- 
ber fluctuation but the much enhanced phase fluctuation 
than the wave function without Gutzwiller projection in 
the underdoped region. We also found AN^AQ much 
greater than 1. 

In this paper, we only used a uniform fugacity fac- 
tor in the d-RVB states. Without including a Jastrow 
factor to obtain a lower energy, this factor is close to 
the derived results of GA. However, including the Jas- 
trow factor produces a much larger fugacity factor and a 
signiflcantly enhanced phase fiuctuation. This indicates 
that the fugacity factor may be quite important in the 
underdoped region. In addition, the fugacity factor could 
have a spatial dependence or a momentum dependence 
as noticed by Anderson'^^ . The spatial dependence could 
produce the stripes as we demonstrated in Refs. [13 and 
l44lthat the Gutzwiller projection operator introduces the 
coupling between charge density, spin density and pair 
fields. The effect of the momentum dependence will be 
left for future study. 
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